Higher‐order functional connectivity analysis of resting‐state functional magnetic resonance imaging data using multivariate cumulants

Abstract Blood‐level oxygenation‐dependent (BOLD) functional magnetic resonance imaging (fMRI) is the most common modality to study functional connectivity in the human brain. Most research to date has focused on connectivity between pairs of brain regions. However, attention has recently turned towards connectivity involving more than two regions, that is, higher‐order connectivity. It is not yet clear how higher‐order connectivity can best be quantified. The measures that are currently in use cannot distinguish between pairwise (i.e., second‐order) and higher‐order connectivity. We show that genuine higher‐order connectivity can be quantified by using multivariate cumulants. We explore the use of multivariate cumulants for quantifying higher‐order connectivity and the performance of block bootstrapping for statistical inference. In particular, we formulate a generative model for fMRI signals exhibiting higher‐order connectivity and use it to assess bias, standard errors, and detection probabilities. Application to resting‐state fMRI data from the Human Connectome Project demonstrates that spontaneous fMRI signals are organized into higher‐order networks that are distinct from second‐order resting‐state networks. Application to a clinical cohort of patients with multiple sclerosis further demonstrates that cumulants can be used to classify disease groups and explain behavioral variability. Hence, we present a novel framework to reliably estimate genuine higher‐order connectivity in fMRI data which can be used for constructing hyperedges, and finally, which can readily be applied to fMRI data from populations with neuropsychiatric disease or cognitive neuroscientific experiments.


A. Asymptotic sampling distributions of coskewness and cokurtosis
Let X ∈ R n be a zero-mean random vector modeling the fMRI signals from n brain regions.The coskewness between regions i, j, k is defined as and its plug-in estimator is obtained by replacing the expectations by sample averages over a random sample of size N : We derive the asymptotic distribution of rc i,j,k .Define the random vector and let θ and Ω be, respectively, the expected value and covariance matrix of Y .Application of the central limit theorem to the sample mean ⟨Y ⟩ shows that √ N (⟨Y ⟩ − θ) ∼ N (0, Ω) for large N .We now apply Cramer's theorem to the function g(x, y, z, u) = x (yzu) 3/2 , which shows that √ N (g(⟨Y ⟩) − g(θ)) ∼ N (0, ∇ g (θ) ′ Ω∇ g (θ)), where ∇ g (θ) is the gradient of g at θ. But and therefore g(θ) = r c i,j,k, .Furthermore, g(⟨Y ⟩) = rc i,j,k so that ) for large N .This shows that rc i,j,k is a consistent estimator of r c i,j,k and is asymptotically normal with variance σ 2 r = ∇ g (θ) ′ Ω∇ g (θ)/N .

B. Asymptotic sampling distribution of the edge connectivity
Let X ∈ R n be a zero-mean random vector and let f 1 (X) and f 2 (X) be monomials in the components of X. Below we suppress the dependence of f 1 and f 2 on X.The edge connectivity has the following form: , for particular choices of f 1 and f 2 , and its plug-in estimator is defined by replacing expectations by sample averages over a random sample of size N : We derive the distribution of r for large N .Define the random vector Y ∈ R 3 by where ′ denotes the transpose, and let θ and Ω be, respectively, the expected vector and covariance matrix of Y .Application of the central limit theorem to the sample mean ⟨Y ⟩ shows that √ N (⟨Y ⟩ − θ) ∼ N (0, Ω) for large N .We now apply Cramer's theorem to the function g(x, y, z) = x/ √ yz, which , where ∇ g (θ) is the gradient of g at θ.But ′ hence g(θ) = r, and furthermore, g(⟨Y ⟩) = r, so that √ N (r − r) ∼ N (0, ∇ g (θ) ′ Ω∇ g (θ)) for large N , where This shows that r is a consistent estimator of r and is asymptotically normal with variance σ 2 r = ∇ g (θ) ′ Ω∇ g (θ)/N .

C. Example fMRI signals with strong higher-order interactions
We selected the three subjects with the largest absolute average coskewness with the primary motor cortices and subsequently selected the fMRI signals from the primary motor cortices and from the target region with the largest absolute coskewness.The signals are shown in Figure 1.We indeed observe coherent extreme deactivations of all three regions are several moments of the scanning session, for example at around 200 seconds in HCP subject 9.

D. Higher-order connectivity maps in single-subject data
To assess if higher-order connectivity can be detected in single-subject data, we used one of the four homologous region-pairs from Section 3.2.2 in the main text, namely the frontal eye-fields.As in Section 3.2.2 of the main text, we tested for statistical significance of third-order connectivity with all cortical regions and for fourth-order connectivity with all homologous region-pairs.The resulting seed-based maps can therefore directly be compared to those obtained on the group-level.Testing was done using the block bootstrap with block-length L = 10 samples by computing 95% confidence intervals, based on 1000 bootstrap realizations.In contrast to the simulated data, the connectivity estimates were not entirely normally distributed, so that a Gaussian approximation to the sampling distributions could not be used to obtain the intervals.The intervals were therefore computed by taking the 2.5% and 97.5% percentiles of the sampling distributions.
The raw seed-based maps are displayed in Figure 1A (coskewness), B (cokurtosis) and C (edge connectivity).Note the large inter-subject variability in the maps.Specifically, the average spatial correlation between the maps from different subjects is 0.08 ± 0.24 (coskewnss), 0.02 ± 0.19 (cokurtosis), and 0.01 ± 0.14 (edge connectivity).The simulations carried out in Section 3.1.2of the main text suggest that this variability reflects statistical uncertainty, rather than variability in the true connectivity maps.This is also suggested by the numerical values of the connectivity estimates (see the colorbars in Figure 1A-C), which are about an order of magnitude larger than those obtained on the group-level (see Figure 6 of the main text).To confirm this, we computed the spatial correlations of the maps with those obtained from a separate scanning session.The subject-averaged spatial correlations were 0.17 ± 0.30 (coskewness), 0.10 ± 0.23 (cokurtosis), and 0.08 ± 0.17 (edge connectivity).
These low values show that the intra-subject variability in the connectivity maps is particularly large.
For comparison, the subject-averaged spatial correlation between the second-order correlations maps obtained from two separate scanning sessions is 0.66 ± 0.12.From this analysis we therefore conclude that the large inter-subject variability in the third-and fourth-order connectivity maps indeed reflects statistical uncertainty.
The thresholded seed-based maps are shown in Figure 1D (coskewness), E (cokurtosis), and F (edge connectivity).The figures show, for each cortical region, the number of subjects (out of a total of 94 subjects) for which the fMRI signal from that region was significantly connected to the two reference regions (left and right frontal eye-fields).The maps resemble the corresponding group-level maps in Figures 7 (coskewness) and 9 (cokurtosis and edge connectivity) of the main text.This shows that the regions that constitute the group-level maps are more likely to be statistically significant in single-subject data than other cortical regions.However, these regions are significant only in 26 (skewness), 19 (kurtosis), and 18 (edge connectivity) subjects.This analysis thus shows that thirdand fourth-order connectivity can be detected in single-subject data, but only in a minority of the subjects.Furthermore, in the above statistical analysis, no correction for multiple corrections was applied.Such a correction requires the approximation of extreme percentiles of the estimators' sampling distributions, namely 2.5/360 = 0.0069% and 100 − 97.5/360 = 99.7292%,which requires an unfeasibly large number of bootstrap realizations.When the sampling distributions of the plug-in estimators are approximated by normal distributions, the statistical threshold is set to α = 0.01, and a Bonferroni-correction is applied, we found no significant third-and fourth-order connectivity in any of the 94 subjects.To check if the results are substantially different when using different regions, we repeated the above analysis using the precunei as seed regions, instead of the frontal eye-fields.The results were comparable with those for the frontal eye-fields.Specifically, the average spatial correlation between the maps from different subjects is 0.03 ± 0.23 (coskewnss), 0.03 ± 0.15 (cokurtosis), and 0.02 ± 0.13 (edge connectivity).We again computed the spatial correlations of the maps with those obtained from a separate scanning session.The subject-averaged spatial correlations were 0.14 ± 0.26 (coskewness), 0.06 ± 0.18 (cokurtosis), and 0.15 ± 0.15 (edge connectivity).For comparison, the subject-averaged spatial correlation between the second-order correlations maps obtained from two separate scanning sessions is 0.67 ± 0.12.Lastly, the maximum number of subjects for which regions were significant are 27 (skewness), 20 (kurtosis), and 18 (edge connectivity).

E. Symmetric random variables have vanishing odd multivariate moments
We show that the multivariate moments of odd order of three random variables with symmetric marginal distributions are zero.In particular, the coskewness of three symmetric random variables is zero.Let X, Y , and Z be zero-mean random variables with joint probability distribution f (x, y, z) and suppose that their marginal distributions are even, i.e.
f (x, y, z)dxdy = f (x, y, −z)dxdy, and the same for x and y and let be the third moment of (x, y, z).This implies that xyf (x, y, z)dxdy = xyf (x, y, −z)dxdy, for all functions g(x, y) and the same for x and y.Substituting (x ′ , y ′ , z ′ ) = (x, y, z) shows that Note that this step uses the oddness of the moment.Subsequently applying Eq. ( 15) to Eq. ( 16) and using Eq. ( 12), i.e. for each of the three variables, shows that ⟨XY Z⟩ = − ⟨XY Z⟩ and hence ⟨XY Z⟩ = 0.The general case is proved in the same way.
F. Finite moments of resting-state fMRI signals Note that the sampling variability of the moments is large for small sample sizes, but stabilizes when the sample size becomes larger, indicating that the corresponding theoretical moments are finite.

Figure 1 :
Figure 1: fMRI signals with strong third-order correlation.fMRI signals from the left (red) and right(blue) primary motor cortex, and a target cortical region (yellow) for three different HCP subjects.For each of the three subjects, the target region was chosen so that its average third-order correlation with the fMRI signals in the primary motor cortices was maximal in absolute value.For better visibility, all fMRI signals were normalized to unit variance.

Figure 2 :
Figure 2: fMRI signals with strong fourth-order correlation.fMRI signals from the right (blue) and left (red) primary motor cortex, together with right (yellow) and left (purple) target regions for three different HCP subjects.For each of the three subjects, the target regions were chosen so that their average fourth-order correlation with the fMRI signals in the primary motor cortices was maximal in absolute value.For better visibility, all fMRI signals were normalized to unit variance.

Figure 1 :
Figure 1: Higher-order connectivity in single-subject data.A. Coskewness maps for each of the 94 HCP subjects.B. Cokurtosis maps for each of the 94 HCP subjects.C. Edge connectivity for each of the 94 HCP subjects.D. Number of subjects (out of the 94 subjects) for which the coskewnesswith the frontal eye fields was significant.E. and F. Same format as D but for the cokurtosis and edge connectivity, respectively.All maps were seeded at the left and right frontal eye fields.For the coskewness maps, the target varied over all cortical regions, and for the cokurtosis and edge connectivity maps, the target varied over all homologous region-pairs.The significance threshold was set to α = 0.05 (uncorrected).

Figure 1 :
Figure 1: Sample moments as a function of sample size.Shown are the first sample moment (i.e.the sample mean), the second-order sample moment (i.e. the sample covariance), the third-order sample moment and the fourth-order sample moment, as functions of sample size, for 10 4 randomly selected quadruples of regions of a single HCP subject.The sample size ranges from 1 to 4800, which corresponds to four scanning sessions (each HCP subject has four scanning sessions of 1200 samples).